function [dy algebraic] = CICR_func(t,y,CONSTANTS)
global nAlgebraic shouldCalculateAlgebraic
% From: http://www.ncbi.nlm.nih.gov/pubmed/15142845
% Title: Computer model of action potential of mouse ventricular myocytes
% Authors: Bondarenko, Szigeti, Bett, Kim, Rasmusson, 2004

% global V_m C_m vol_myo vol_JSR  vol_NSR vol_ss A_cap Ca_o R T F CMDN_tot CSQN_tot Km_CMDN Km_CSQN ...
%     v1 v2 v3 tau_tr tau_xfer K_m_up  i_CaL_max k_plus_a k_minus_a k_plus_b ...
%     k_minus_b k_plus_c k_minus_c m n E_CaL G_CaL G_CaB K_pcb k_pc_max k_pc_half G_Ca_total stimVoltage...
%     Na_i Na_o k_NaCa K_mNa K_mCa k_sat eta i_pCa_max Km_pCa


% -------------------------CONSTANTS----------------------------
V_m  = CONSTANTS(1);
C_m = CONSTANTS(2);
vol_myo  = CONSTANTS(3);
vol_JSR  = CONSTANTS(4);
vol_NSR = CONSTANTS(5);
vol_ss  = CONSTANTS(6);
A_cap  = CONSTANTS(7);
Ca_o  = CONSTANTS(8);
R  = CONSTANTS(9);
T = CONSTANTS(10);
F = CONSTANTS(11);
CMDN_tot  = CONSTANTS(12);
CSQN_tot = CONSTANTS(13);
Km_CMDN  = CONSTANTS(14);
Km_CSQN  = CONSTANTS(15);
v1  = CONSTANTS(16);
v2  = CONSTANTS(17);
v3  = CONSTANTS(18);
tau_tr  = CONSTANTS(19);
tau_xfer = CONSTANTS(20);
K_m_up = CONSTANTS(21);
i_CaL_max = CONSTANTS(22);
k_plus_a = CONSTANTS(23);
k_minus_a = CONSTANTS(24);
k_plus_b = CONSTANTS(25);
k_minus_b  = CONSTANTS(26);
k_plus_c = CONSTANTS(27);
k_minus_c = CONSTANTS(28);
m  = CONSTANTS(29);
n = CONSTANTS(30);
E_CaL = CONSTANTS(31);

G_CaB = CONSTANTS(32);
K_pcb = CONSTANTS(33);
k_pc_max = CONSTANTS(34);
k_pc_half = CONSTANTS(35);
G_CaT = CONSTANTS(36);
Na_i  = CONSTANTS(37);
Na_o = CONSTANTS(38);
k_NaCa = CONSTANTS(39);
K_mNa = CONSTANTS(40);
K_mCa = CONSTANTS(41);
k_sat = CONSTANTS(42);
eta  = CONSTANTS(43);
i_pCa_max = CONSTANTS(44);
Km_pCa = CONSTANTS(45);
G_CaL = CONSTANTS(46); 
% RyR Constants
koCa = CONSTANTS(47);
kom  = CONSTANTS(48); 
kiCa  = CONSTANTS(49);
kim  = CONSTANTS(50);
EC50_SR  = CONSTANTS(51);
ks  = CONSTANTS(52);
MaxSR  = CONSTANTS(53);
MinSR  = CONSTANTS(54);
HSR = CONSTANTS(55);

% ---------------------------------------------------------------------


global TMAX

V_m(t>TMAX/2) = -30;  
V_m(t<TMAX/2) = -60; 
% if t > TMAX/2
%     V_m  = -30;
% else
%     V_m =  -60;
% end


O = y( 1,:);
C2= y(2,:);
C3= y(3, :);
C4= y( 4, :);
%     --------RyR------- debug---------------
% P_C2= y( 5, :);
% P_O1= y( 6,:);
% P_O2= y( 7,:);


    ryr_R = y(5);
    ryr_O = y(6);
    ryr_I = y(7);
    
%     -----------------------

I1= y( 8,:);
I2= y( 9,:);
I3= y( 10,:);
Ca_JSR= y( 11,:);
Ca_NSR= y( 12,:);
Ca_ss= y( 13,:);
Ca_i= y( 14,:);
P_RyR = y( 15,:);
b = y(16,:);
g = y( 17,:);
% ---RyR----- debug-----
ryr_RI = y(18);

% -----------------------------------------------------------------------


% Component: Ryanodine Receptor

% P_C1 = 1.00000 - (P_C2+P_O1+P_O2);
% P_O1_prime = ( k_plus_a .* (Ca_ss .^ n) .* P_C1+ k_minus_b .* P_O2+ k_minus_c .* P_C2) - ( k_minus_a .* P_O1+ k_plus_b .* (Ca_ss .^ m) .* P_O1+ k_plus_c .* P_O1);
% P_O2_prime =   k_plus_b.*(Ca_ss .^ m) .* P_O1 -  k_minus_b.*P_O2;
% P_C2_prime =   k_plus_c.*P_O1 -  k_minus_c.*P_C2;

CajSR = Ca_JSR/1e3; % convert mM to uM
Casub = Ca_ss/1e3; % convert mM to uM

jSRCarel = ks .* ryr_O .* (CajSR - Casub) ;
kCaSR = MaxSR - (MaxSR - MinSR)./ (1 + (EC50_SR./CajSR).^HSR);
koSRCa = koCa./kCaSR;
kiSRCa = kiCa.*kCaSR;
d_ryr_R_prime = (kim.*ryr_RI - kiSRCa.*Casub.*ryr_R) - (koSRCa .* Casub.^2.*ryr_R - kom.*ryr_O) ;
d_ryr_O_prime =(koSRCa.*Casub.^2 .*ryr_R - kom.*ryr_O) - (kiSRCa .* Casub.*ryr_O - kim.*ryr_I) ;
d_ryr_I_prime = (kiSRCa.* Casub.*ryr_O - kim.*ryr_I) - (kom.*ryr_I - koSRCa.*Casub.^2 .* ryr_RI) ;
d_ryr_RI_prime = (kom.*ryr_I - koSRCa.*Casub.^2.*ryr_RI)-(kim.*ryr_RI - kiSRCa.*Casub.*ryr_R);

% Component: L-type Calcium Current

alpha = ( 0.400000.*(exp(((V_m+12.0000)./10.0000))).*((1.00000+ 0.700000.*(exp(( - ((V_m+40.0000) .^ 2.00000)./10.0000)))) -  0.750000.*(exp(( - ((V_m+20.0000) .^ 2.00000)./400.000)))))./(1.00000+ 0.120000.*(exp(((V_m+12.0000)./10.0000))));
beta =  0.0500000.*(exp(( - (V_m+12.0000)./13.0000)));
gamma = ( k_pc_max.*Ca_ss)./(k_pc_half+Ca_ss);
K_pcf =  13.0000.*(1.00000 - (exp(( - ((V_m+14.5000) .^ 2.00000)./100.000))));

C1 = 1.00000 - (O+C2+C3+C4+I1+I2+I3);
C2_prime = ( 4.00000.*alpha.*C1+ 2.00000.*beta.*C3) - ( beta.*C2+ 3.00000.*alpha.*C2);
C3_prime = ( 3.00000.*alpha.*C2+ 3.00000.*beta.*C4) - ( 2.00000.*beta.*C3+ 2.00000.*alpha.*C3);
I1_prime = ( gamma.*O+ 0.00100000.*( alpha.*I3 -  K_pcf.*I1)+ 0.0100000.*( alpha.*gamma.*C4 -  4.00000.*beta.*K_pcf.*I1)) -  K_pcb.*I1;
I2_prime = ( 0.00100000.*( K_pcf.*O -  alpha.*I2)+ K_pcb.*I3+ 0.00200000.*( K_pcf.*C4 -  4.00000.*beta.*I2)) -  gamma.*I2;
I3_prime = ( 0.00100000.*( K_pcf.*I1 -  alpha.*I3)+ gamma.*I2+ 1.00000.*gamma.*K_pcf.*C4) - ( 4.00000.*beta.*K_pcb.*I3+ K_pcb.*I3);
O_prime = ( alpha.*C4+ K_pcb.*I1+ 0.00100000.*( alpha.*I2 -  K_pcf.*O)) - ( 4.00000.*beta.*O+ gamma.*O);
C4_prime = ( 2.00000.*alpha.*C3+ 4.00000.*beta.*O+ 0.0100000.*( 4.00000.*K_pcb.*beta.*I1 -  alpha.*gamma.*C4)+ 0.00200000.*( 4.00000.*beta.*I2 -  K_pcf.*C4)+ 4.00000.*beta.*K_pcb.*I3) - ( 3.00000.*beta.*C4+ alpha.*C4+ 1.00000.*gamma.*K_pcf.*C4);


I_CaL =  G_CaL.*O.*(V_m - E_CaL); % picoA
B_ss = (1.00000+( CMDN_tot.*Km_CMDN)./((Km_CMDN+Ca_ss) .^ 2.00000)) .^ ( - 1.00000);

% NCX current
i_NaCa =(( (( (( k_NaCa.*1)./((K_mNa .^ 3)+(Na_o .^ 3))).*1)./(K_mCa+Ca_o)).*1)./(1+ k_sat.*(exp((( (eta - 1).*V_m.*F)./( R.*T)))))).*( (exp((( eta.*V_m.*F)./( R.*T)))).*(Na_i .^ 3).*Ca_o -  (exp((( (eta - 1).*V_m.*F)./( R.*T)))).*(Na_o .^ 3).*Ca_i);

% PMCA current
i_pCa = ( i_pCa_max.*(Ca_i .^ 2))./((Km_pCa.^ 2)+(Ca_i .^ 2));


% Background Ca current
E_CaN =  (( R.*T)./( 2.*F)).*(log((Ca_o./Ca_i)));
i_CaB = G_CaB.*(V_m - E_CaN); % (microA_per_microF)

% Voltage dependent T-type calcium current I_CaT
% The T-type model has been adopted from the paper, "IONIC MECHANISMS FOR
% ELECTRICAL HETEROGENEITY BETWEEN RABBIT PURKINJE AND VENTRICULAR CELLS,
% Oleg V. Aslanidi et al.

g_inf = 1.0./(1+ exp((V_m + 60.0)./6.6));
alpha_g = 0.015.*exp(-(V_m+71.7)./83.33);
beta_g = 0.015.*exp((V_m+71.7)./15.38);
tau_g = 1.0./(alpha_g+ beta_g);

g_prime = (g_inf -g) ./ tau_g;
alpha_b = 1.068.*exp((V_m + 16.3)./30.0);
beta_b =  1.068.*exp(-(V_m + 16.3)./30.0);

b_inf = 1.0./(1+ exp(-(V_m + 28.0)./6.1));
tau_b = 1.0./(alpha_b+ beta_b);

b_prime = (b_inf - b )./tau_b;

I_CaT = G_CaT .*b.*g.*(V_m - 50); % micro A

% Calcium Current of both T-type and L-type Ca Channels
I_Ca_channels= I_CaT + I_CaL;

% Component: Calcium Fluxes
% J_rel =  v1.*(P_O1+P_O2).*(Ca_JSR - Ca_ss).*P_RyR; % (micromolar_per_millisecond)
J_rel = jSRCarel * 1e3; % debug

J_tr = (Ca_NSR - Ca_JSR)./tau_tr; %  (micromolar_per_millisecond)
J_leak =  v2.*(Ca_NSR - Ca_i); %  (micromolar_per_millisecond)
J_up = ( v3.*(Ca_i .^ 2))./((K_m_up .^ 2)+(Ca_i .^ 2)); %  (micromolar_per_millisecond) -
J_xfer = (Ca_ss - Ca_i)./tau_xfer; %  (micromolar_per_millisecond)
P_RyR_prime =  - 0.04.*P_RyR -  (( 0.1.*I_Ca_channels)./i_CaL_max).*(exp(( - ((V_m - 5) .^ 2)./648))); %  (micromolar_per_millisecond)

% Component: Calcium Concentration
B_i = (1+( CMDN_tot.*Km_CMDN)./((Km_CMDN+Ca_i) .^ 2)) .^ ( - 1); % Calmodulin buffer

% Ca_i_prime =  B_i.*((J_leak+J_xfer) - (J_up+( ((i_CaB) ).*A_cap.*C_m)./( 2.*vol_myo.*F))); % I deleted the 'J_trpn', NCX and PMCA currents terms as we are not interested in it yet
Ca_i_prime =  B_i.*((J_leak+J_xfer) - (J_up+( ((i_CaB+i_pCa) -  2.*i_NaCa).*A_cap.*C_m)./( 2.*vol_myo.*F)));
Ca_ss_prime =  B_ss.*(( J_rel.*vol_JSR)./vol_ss - (( J_xfer.*vol_myo)./vol_ss+( I_Ca_channels.*A_cap.*C_m)./( 2.*vol_ss.*F)));
B_JSR = (1+( CSQN_tot.*Km_CSQN)./((Km_CSQN+Ca_JSR) .^ 2)) .^ ( - 1);
Ca_JSR_prime =  B_JSR.*(J_tr - J_rel);
Ca_NSR_prime = ( (J_up - J_leak).*vol_myo)./vol_NSR - ( J_tr.*vol_JSR)./vol_NSR;

%--------------------algebraic variables------------------------
if shouldCalculateAlgebraic
    algebraic = zeros(length(t),8);
    algebraic(:,1) = I_CaL';
    algebraic(:,2) = i_NaCa';
    algebraic(:,3) = i_pCa';
    algebraic(:,4) = i_CaB';
    algebraic(:,5) = I_CaT';
    algebraic(:,6) = J_rel';
    algebraic(:,7) = J_up';
    algebraic(:,8) = V_m';
end
%------------------------rates-----------------------------------------
dy = zeros(length(y),1);
if  ~shouldCalculateAlgebraic
    dy(1) = O_prime;
    dy(2) = C2_prime;
    dy(3) = C3_prime;
    dy(4) = C4_prime;
%     --------RyR------- debug---------------
%     dy(5) = P_C2_prime;
%     dy(6) = P_O1_prime;
%     dy(7) = P_O2_prime;
  dy(5) = d_ryr_R_prime;
  dy(6) = d_ryr_O_prime;
  dy(7) = d_ryr_I_prime;
%   ---------------------------------------------
    dy(8) = I1_prime;
    dy(9) = I2_prime;
    dy(10) = I3_prime;
    dy(11) = Ca_JSR_prime;
    dy(12) = Ca_NSR_prime;
    dy(13) = Ca_ss_prime;
    dy(14) = Ca_i_prime;
    dy(15) = P_RyR_prime;
    dy(16) = b_prime;
    dy(17) = g_prime;
    %     --------RyR------- debug---------------
 dy(18) = d_ryr_RI_prime;
end
end